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We present a general multi-component density functional theory in which electrons and nuclei are 
treated completely quantum mechanically, without the use of a Born-Oppenheimer approximation. 
The two fundamental quantities in terms of which our theory is formulated are the nuclear N-body 
density and the electron density expressed in coordinates referring to the nuclear framework. For 
these two densities coupled Kohn-Sham equations are derived and the electron-nuclear correlation 
functional is analyzed in detail. The formalism is tested on the hydrogen molecule H2 and its 
positive ion H2 using several approximations for the electron-nuclear correlation functional. 
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I. INTRODUCTION 



Density functional theory (DFT) is among the most 
succesful approaches to calculate the electronic structure 
of atoms, molecules and solids. In its original form 
DFT always invokes the Born-Oppenheimer approxima- 
tion: One is supposed to calculate the electron density 
p{r) which is in 1-1 correspondence to the static poten- 
tial potential of fixed nuclei. In a recent Letter Q we 
introduced a multicomponent density-functional theory 
(MCDFT) for the complete quantum treatment of many- 
particle systems consisting of electrons and nuclei. With 
this theory it is possible to decribe from first principles 
physical phenomena that depend on a strong coupling 
between electronic and nuclear motion. MCDFT thereby 
extends the widely applied density functional formalism 
for purely electronic properties, opening up a new field 
of applications, such as the first-principles calculation of 
electron-phonon coupling in solids 0| which is a key in- 
gredient in the description of superconductivity 
and polaronic motion 9, 10]. The quantum treatment 
of the nuclear motion in molecules or solids is essential 
in situations that from a Born-Oppenheimer (BO) view- 
point must be described by a superposition of different 
BO structures. This is, for instance, the case in flo ppy 
molecules so-called switchable molecules |l2j 

which are in a superposition of an open and a closed state 
after a laser excitation. Apart from treating such various 
phenomena the MCDFT presented here also paves the 
way for future time-dependent extensions of the theory 
which would enable one to calculate the coupled elec- 
tronic and nuclear dynamics of many-particle systems, 
within linear response and beyond. Indeed, some prelim- 
inary steps towards the description of the coupled ion- 
ization and dissociation dynamics of molecules in strong 
laser fields have already been taken |0, 0, ^ . 
The purpose of the present work is twofold. First, we 



want to give an extended and detailed description of the 
theory that was briefly described in our Letter. Second, 
we want to investigate in detail some new approximate 
density functionals for the electron-nuclear correlation 
and see how they perform. To do this the formalism 
is tested on the hydrogen molecule and its positive ion. 
The paper is organized as follows. In section II we first 
introduce the basic formalism and discuss the Hohcnbcrg- 
Kohn theorem and the Kohn-Sham equations in a multi- 
component theory in which the electron density is defined 
with respect to a coordinate frame attached to the nu- 
clear framework and in which the diagonal of the nuclear 
density matrix appears as a new variable. In section III 
we perform an analysis of the several energy functionals 
and of the resulting potentials in the Kohn-Sham equa- 
tions. Furthermore, the connections between the effec- 
tive potential of the nuclear Kohn-Sham equation and the 
Born-Oppenheimer energy surface is analyzed. In section 
IV we apply our formalism and test several approximate 
forms for the electron-nuclear correlation functional for 
the case of the hydrogen molecule and its positive ion. 
Finally in section V we present our conclusions. 



II. BASIC FORMALISM 

A. Discussion of the Hamiltonian 

We consider a system composed of Ng electrons with 
coordinates {vj} = r and Nn nuclei with masses 
Mi...Mn„, charges Zi...Zn^, and coordinates denoted by 
{Rq} = R. By convention, the subscript "e" and "n" re- 
fer to electrons and nuclei, respectively, and atomic units 
are employed throughout this work. In non-relativistic 
quantum mechanics, the system is decribed by the Hamil- 
tonian 

H - f„(R) + Ty„„(R)-HC/oxt,„(R) 



2 



where 



where 



(1) 

(2) 
(3) 



denote the kinetic-energy operators of the nuclei and elec- 
trons, respectively, and 



We, 



We, 



ZnZ, 



2 j-^^ |R-a - R-/3I 



i,i=\ 

EE' 



(4) 



(5) 



(6) 



represent the interparticle Coulomb interactions. We em- 
phasize that no BO approximation has been assumed 
in the Hamiltonian of Eq.l^ provides a quantum 
mechanical description of all, i.e., electronic and nuclear 
degrees of freedom. In contrast to the standard approach 
using the BO approximation, the interactions between 
electrons and nuclei are therefore treated within Wen, 
Eq. ijHJl, and do not contribute to the external potentials. 
Truly external potentials representing, e.g., a voltage ap- 
plied to the system, are contained in 



N,^ 



C^oxt,n(R-Q) 



a=l 



(7) 

(8) 



Defining electronic and nuclear single-particle densities 
conjugated to the true external potential iQ, a MCDFT 
formalism can readily be formulated on the basis of the 
above Hamiltonian However, as discussed in 

such a MCDFT is not useful in practice because the 
single-particle densities necessarily reflect the symmetry 
of the true external potentials and are therefore not char- 
acteristic of the internal properties of the system. In par- 
ticular, for all isolated systems where the external poten- 
tials (TJ vanish, these densities, as a consequence of the 
translational invariance of the respective Hamiltonian, 
are constant. 

A suitable MCDFT is obtained by defining the den- 
sities with respect to internal coordinates of the system 
To this end, new electronic coordinates are intro- 
duced according to 



r'j = 7^(a, /3, 7) (r^ 



R, 



CMN 



l...iV, 



(9) 



R 



■CMN 



(10) 



a=l 



denotes the center of mass (CM) of the nuclei, the total 
nuclear mass is given by 



a=l 



(11) 



and TZ is the three-dimensional orthogonal matrix rep- 
resenting the Euler rotations The Euler angles 
(a, /9, 7) are functions of the nuclear coordinates {R} 
and specify the orientation of the body-fixed coordi- 
nate frame. They can be determined in various ways. 
One way to define them is by requiring the inertial 
tensor of the nuclei to be diagonal in the body-fixed 
frame. The conditions that the off-diagonal elements of 
the inertia tensor are zero in terms of the rotated co- 
ordinates 7^(Rq — Rcj\/A') then give three determining 
equations for the three Euler angles in terms of the nu- 
clear coordinates {R}. This way of choosing the Eu- 
ler angles is commonly used within the field of nuclear 
physics 0, 0, li^l but is, of course, not unique. A 
common alternative way to determine the orientation of 
the body-fixed s^tem is provided by the so-called Eckart 
conditions [H I2I IH 111 IH IH [23 which are suitable 
to describe small vibrations in molecules and phonons in 
solids 3. A general and very elegant discussion on the 
various ways the body-fixed frame can be chosen is given 
in reference {28] . In this work we will not make a specific 
choice as our derivations are independent of such a choice. 
The most important point is that by virtue of Eq. 
the electronic coordinates are defined with respect to a 
coordinate frame that is attached to the nuclear frame- 
work and rotates as the nuclear framework rotates. In 
fact, this transformation comprises two transformations: 
A first one transforming the space-fixed inertial coordi- 
nates into CM-fixed relative coordinates, and a second 
one transforming the CM-fixed relative coordinates into 
body- fixed internal coordinates. 

The nuclear coordinates themselves are not trans- 
formed any further at this point, i.e., 



Rl — Rn 



l...N„. 



(12) 



Of course, introducing internal nuclear coordinates is also 
desirable. However, the choice of such coordinates de- 
pends strongly on the specific system to be described: If 
near-equilibrium situations in systems with well-defined 
geometries are considered, normal or - for a solid - 
phonon coordinates are most appropriate, whereas frag- 
mentation processes of molecules are better described in 
terms of Jacobi coordinates Therefore, keeping a 

high degree of fiexibility, the nuclear coordinates are left 
unchanged for the time being and are only transformed 
to internal coordinates prior to actual applications in the 
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final equations that we will derive. Another reason for 
not introducing any internal nuclear coordinates at this 
point, is to retain simple forms of the equations. In a 
transformation to internal nuclear coordinates typically 
the nuclear center-of-mass and the Eulcr angles are taken 
as new variables as well as 3A^„ — 6 internal or shape co- 
ordinates Qi I27I [^. These internal coordinates, 
however, do not have a simple relation to the original Nn 
nuclear coordinates and will therefore lead to a compli- 
cated form of the Hamiltonian in the new coordinates. 
We will therefore delay the use of such transformations 
until we have derived the final equations. 

As a result of the coordinate changes of Eq.© , the 
Hamiltonian |^ transforms into 

+ fe(l') + W^ee(l') + TmPC®,I') 

+ I^en(a,l') + C>cxt,e(E,l')- (13) 

Since we have transformed to a noninertial coordinate 
frame mass-polarization and Coriolis (MPC) terms 




(14) 

appear. Obviously, Tmpc is not symmetric in the elec- 
tronic and nuclear coordinates. However, this was not 
expected since only the electrons refer to a noninertial co- 
ordinate frame, whereas the nuclei are still defined with 
respect to the inertial frame. Therefore, all MPC terms 
arise solely from the electronic coordinates, representing 
ficticious forces due to the electronic motion in noniner- 
tial systems (for a detailed form of these terms within the 
current coordinate transformation see (H) . The kinetic- 
energy operators Tg and T„, the electron-electron and 
nuclear- nuclear interactions, as well as the true exter- 
nal potential Uext,n acting on the nuclei are formally 
unchanged in Eq. H13|l and therefore given by Eqs. Q 
and Q with the new coordinates replacing the old ones, 
whereas the electron-nuclear interaction now reads 

" |r;.-7^(a,/3,7)(R.-RcAfA^)|■ ^^^^ 

The quantity 

R^=7^(a,/3,7)(R„-Rc•A^iv) (16) 

that appears in Ea. (|15|l is a so-called shape coordinate 4, 
|2^ . i.e. it is invariant under rotations and translations 
of the nuclear framework 

R::(OE + a)=R:(E) (17) 



where O is an arbitrary rotation matrix and a an ar- 
bitrary translation vector. The invariance property de- 
scribed in Eq. H17() is simply a consequence of the fact that 
the Eulcr angles are defined by giving the vectors R" 
certain values, independent of where the nuclear center- 
of-mass was situated in the laboratory frame or how the 
nuclear framework was orientated. This is, of course, 
precisely the purpose of introducing a body-fixed frame. 
For this reason the potential in Ea. (|15|l that the elec- 
trons in the body-fixed frame experience from the nuclei 
is invariant under rotations or translations of the nuclear 
framework. 

As a further result of the coordinate transformation 
0, the true external potential acting on the electrons 
now not only depends on the electronic coordinates, but 
also on all the nuclear coordinates: 

We 

J7ext,e(E,l') = ^Mcxt,e(7^-'r;■ -|- Rc M n) ■ (18) 
J = l 

In the chosen coordinate system the electron-nuclear in- 
teraction H15(l and the external potential (|18|l remain one- 
body operators with respect to the electronic degrees of 
freedom but represent complicated A^„-body interactions 
with respect to the nuclei. We finally discuss some gen- 
eral aspects of our coordinate transformation. If we con- 
sider the symmetry properties of our original Hamilto- 
nian of Eq.JQ) in the absence of external potentials, we 
see that it is invariant under simultaneous translations 
and rotations of all particles, i.e. of both electrons and 
nuclei. This is not true anymore for our transformed 
Hamiltonian. Since we transformed the electronic co- 
ordinates to a body-fixed frame we find that in the ab- 
sence of external potentials the transformed Hamiltonian 
of Ea. H13|l is invariant under translations and rotations 
of nuclear coordinates only. The corresponding ground- 
state wavefunction, if it is nondegenerate, will have the 
same invariance. 

Let us next consider the permutational symmetry. The 
ground state wavefunction of the original Hamiltonian of 
Eq.(^ is antisymmetric under the interchange of elec- 
tronic space-spin coordinates and symmetric or antisym- 
metric under interchange of nuclear space-spin coordi- 
nates of nuclei of the same type, depending on whether 
they are bosons or fermions. The ground state wavefunc- 
tion of the transformed Hamiltonian of Eq. H13|) will also 
be antisymmetric with respect to the interchange of elec- 
tronic space-spin coordinates. However, the symmetry 
properties with respect to the interchange of the nuclear 
space-spin coordinates depend on the conditions that we 
choose to determine the Euler angles. If we choose a de- 
termining constraint for the Euler angles that is symmet- 
ric in the interchange of particles of the same type, then 
the transformed wavefunction will retain the permuta- 
tional symmetry properties of the original wavefunction. 
This is, for instance, the case if we determine the Eu- 
ler angles by the requirement that the nuclear inertia 
tensor be diagonal. However, if we choose a nonsym- 
metric constraint, such as the Eckart conditions, then 
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the transformed wavefunction will have more complicated 
transformation properties under the interchange of nu- 
clear spin-space coordinates since the interchange of two 
nuclear coordinates will then also change the Euler an- 
gles (a detailed account on this topic is given in Ref . ) ■ 
This can lead to practical complications but will not af- 
fect our general formalism. 

We finally note that the coordinate transformation we 
presented here did not aim at a separation of the con- 
stants of motion of the system (even for the case of iso- 
lated systems). In contrast, the transformation @ was 
chosen such that the new electronic coordinates reflect 
the internal symmetry of the system. We thus arrive at 
a Hamiltonian which naturally lends itself as a starting 
point for the formulation of a MCDFT, as will be shown 
in the subsequent sections. 



B. Definition of the Densities 

As a first step towards the formulation of a density 
functional theory, one has to define the densities which 
will serve as the fundamental variables of the theory. Al- 
though this seems to be rather straightforward and is 
normally not discussed at length, a careful definition of 
the densities is of crucial importance in the current con- 
text. 

As already mentioned above, it is not useful to define 
electronic and nuclear single-particle densities in terms 
of the inertial coordinates r and R, since such densi- 
ties necessarily reflect the symmetry of the corresponding 
true external potentials, e.g., Galilean symmetry for van- 
ishing external potentials. Therefore, such single-particle 
densities are not characteristic for the internal properties 
of the system under consideration. 

We proceed with the definition of a suitable set of den- 
sities, which should fulfill the following requirements: 

• They should be characteristic for the internal prop- 
erties of the system; in particular, they should be 
meaningful in the limit of vanishing external poten- 
tials. 

• The basic electronic variable should be a single- 
particle quantity. 

• The treatment of the nuclear degrees of freedom 
should allow for appropriate descriptions of situa- 
tions as different as near-equilibrium properties of 
solids and fragmentation processes of molecules. 

A set of densities which meets these requirements is given 
by 

r(E) = ^ / d^-r' |*(E|,i'^)|' (19) 
p(r') = f d^"R f d^^-h' |*(Ei,l'^)|b0) 



where ^(Rs, r'a) corresponds to the ground state of 
Hamiltonian (|13(l and where s and a denote the nuclear 
and electronic spin coordinates. These densities are de- 
fined with respect to the transformed coordinates {R, r'}. 
In particular, the electronic single-particle density p(r') 
refers to the body-fixed molecular frame. In terms these 
coordinates, the quantity H2U|I represents a conditional 
density, which is characteristic for the internal proper- 
ties of the system. It is proportional to the probability 
density of finding an electron at position r' as measured 
from the nuclear center-of-mass, given a certain orienta- 
tion of the nuclear framework. Therefore the electronic 
density calculated through H2U|) reflects the internal sym- 
metries of the system, e.g., the cylindrical symmetry of a 
diatomic molecule, instead of the Galilean symmetry of 
the underlying space. The nuclear degrees of freedom, on 
the other hand, are described using the diagonal of the 
nuclear density matrix, Eq. (|19|l . In the absence of exter- 
nal potentials this quantity will have the transformation 
property 

r(OE + a) = r(E) (21) 

where O is a rotation and a a translation vector. Its per- 
mutational properties will depend on the choice of the 
body fixed frame as discussed in the previous section. 
The quantity r(R) allows us to set up a general as well 
as flexible formalism, which will be applicable to a large 
variety of situations. In an actual application, one may 
at a later stage further contract this quantity to obtain 
reduced density matrices or, depending on the physical 
situation, introduce more suitable internal nuclear coor- 
dinates which could not be done if single-particle quan- 
tities had already been introduced at this point. 

C. The Hohenberg-Kohn Theorem for 
Multicomponent Systems 

In this section, we discuss the extension of the 
Hohenberg-Kohn theorem to multicomponent s :^ste ms. 
In contrast to prior formulations of the MCDFT V 
.31., _3^J , this analysis will employ the densities H19fl 
and as fundamental variables. Correspondingly, the 
starting point of the following analysis is the Hamiltonian 
(|13|1 . In order to formulate a Hohenberg-Kohn- (HK)- 
type statement, the Hamiltonian (|13|l is generalized to 

H = f + W + U + V, (22) 

where 

f = f„ (E) + t (l') + TmpC [R, r') (23) 
denotes the total kinetic-energy operator and 

W ^Wee{i:)+Wen{R,l!) (24) 

contains the electron-electron and the electron-nuclear 
interaction. Furthermore, auxiliary 'external' potentials 
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conjugated to the densities (|19|l and H20() . 

y = K(E) + K(i'), (25) 

have been added to the Hamiltonian. We note that, in 
the transformed coordinates, Vn actually acts as an iV„- 
body operator with respect to the nuclear coordinates. 



(26) 



and particularly contains the internuclear repulsion 
WnnCR), while Ve is a one-body operator with respect 
to the (body- fixed) electronic coordinates: 



ve{r'). 



(27) 



The 'true' external potentials, on the other hand, are 
subsumed in 



C/=C/ext,«(R) + ?/cxt,e(R,r'). 



(28) 



Note that the nuclear potential Ucxt.n has the same struc- 
ture as Vn, whilst the electronic potential Ucxt.e acts sim- 
ilar to the electron-nuclear interaction in the transformed 
coordinate system. 

The Hamiltonian H22|) and the above defined densities 
(|19|l and H20() now provide a suitable basis for the formu- 
lation of the multicomponent Hohenberg-Kohn (MCHK) 
theorem. It can be summarized by the following state- 
ments: 

1. Uniqueness: 

The set of ground-state densities {F, p} uniquely 
determines the ground-state wavefunction, 'J = 
\l/[r_p], as well as the potentials, {Vn = 
V^[r_p], V^e = ye[r_p]}. As a consequence, any ob- 
servable of the static many-body system is a func- 
tional of the set of ground-state densities {F, p}. 

2. MCHK variational principle: 
The total-energy functional 



E[T^p] := (vl/[F,p]|i/|v|/[F,p]) 



(29) 



is equal to the exact ground-state energy Eq if the 
exact densities Fq and po corresponding to fixed ex- 
ternal potentials Vnfi and Ve,o sue inserted into the 
functional. For all other densities, the inequality 



Eo < E[T^p] 



(30) 



holds true. 



This MCHK theorem can be proven by using both 
the reductio ad absurdum and the constrained search ap- 
proach, familiar from standard DFT [s^l- In the fol- 
lowing, a generalization of the latter to multi-component 
(MC) systems will be presented. We start out by defining 
the functional: 



min (*|T 



VF + [/|*), 



(31) 



i.e., we search for the minimum of {^\T + W + C/|\E') 
using all (properly normalized and symmetrized) wave 
functions yielding a given set of densities {F, p}. It must 
be noted that all the wave functions that we use in the 
constrained search procedure are now also required to 
have the correct symmetry properties respect to inter- 
change of nuclear space-spin coordinates of nuclei of the 
same type. As we discussed before these symmetry prop- 
erties depend on the way we define the body-fixed frame. 
For instance, if we define the body-fixed frame by a di- 
agonalization of the nuclear inertia tensor then the con- 
strained search must be carried out over all wavefunctions 
that are antisymmetric in the electronic spin-space coor- 
dinates and symmetric or anti-symmetric with respect 
to the interchange of nuclear spin-space coordinates, de- 
pending on whether the nuclei are bosons or fermions. If 
we denote the minimizing state (assuming it exists _59J) 
by \I>™'"[F,p], we reafize that 

F[T,p] = (^'"""[F.p]|f + iy + C/|^'"""[F./?]) (32) 

is - by construction - a functional of the densities. We 
note that, in contrast to usual DFT, the functional F is 
not universal since it still depends on the external poten- 
tials U which, as a result of our coordinate transforma- 
tion, are functions of both R and r' as was discussed in 
connection with Ea. p8|l . 

Using Ea. (|32|l . the total-energy functional is given by 

E[r^p]^F[r^p]+ J d^"RF(E)K(E) + j dv p[v)v,[v). 

(33) 

The variational principle (|30|l can now be proven by em- 
ploying the Raylcigh-Ritz variational principle: 



Eq = imnmm^) 



(34) 



Following the constrained-search procedure [2S| of ordi- 
nary DFT, the minimum in 134|) is split into two consec- 
utive steps 



^0 



mm I min {■^\H\'^) 
r.p V'l'^r.p 



= mm 
r,P 



F[F,p] + y"d^"RF(E)14(S) 



dv p(r)we(r) 



= min£^[F pi, 
r,p 



(35) 



where the external potentials Vn and v,, are held 
fixed during the minimization (For notational simplic- 
ity, the primes indicating the transformed coordinates are 
dropped from now on. By convention, all electronic coor- 
dinates are understood to refer to the body- fixed frame). 
In the second step, we have exploited the fact that all 
wave functions which lead to the same densities also yield 
the same external energy. By virtue of the Rayleigh- 
Ritz variational principle, the minimizing densities are 
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the ground-state densities Tq and po- Furthermore, any 
other set of densities will lead to an energy above the 
true ground-state energy if inserted in the total-energy 
functional H33|l . This completes the proof of statement 2. 

In order to prove the first statement, we reformulate 
the variational principle (|35|l according to 

d[F[T^p]+ J d^"RT{-^)Vnm + y"drp(rK(r)} =0. 

(36) 

Since the variations can be done independently, Eq. H^fijl 
is equivalent to 



Sp{r) 



VrM) 



(37) 
(38) 



If the exact densities {ro,po} are inserted, the Euler 
equations IpTzjl and are satisfied for the true external 
potentials. If, on the other hand, an arbitrary set of den- 
sities {r, p} is inserted, Eqs. (|37() and (|38|l define - assum- 
ing the functional derivatives exist - a set of potentials, 
which reproduce {F, p} as ground-state densities. There- 
fore, the set of densities {F,p} uniquely determines the 
external potentials {Vn,Ve} and thus the ground-state 
wavefunction ^ = vl/"""[r_p]. 

Before concluding, a number of remarks are added: 

• As usual, the potentials are uniquely determined 
up to an arbitrary additive constant, and non- 
degeneracy of the ground state has been assumed. 

• Similar to purely electronic DFT, the functional 
F[T^p] is defined via Eq. ^ for all {F,p}- 
representable densities, i.e. for all densities ob- 
tained according to Eas. ljiyil and (|20|l from a many- 
body wave function with the right permutational 
symmetries. The potentials {Vn,Ve} are defined 
for all densities, for which the functional derivatives 
in Eqs. (|37|l and (|38|l exist, i.e., for all interacting 
{Vn,Ve}-representable densities. 

• If vanishing external potentials (|18|l are considered, 
the analysis reduces to the one given in |^. 



D. The Kohn-Sham Scheme for Multicomponent 
Systems 

As usual, the HK theorem does not depend on the 
specific form of the particle-particle interaction. In par- 
ticular, it can be applied to an auxiliary system which 
is characterized by = 0, i.e., the system consists of 
noninteracting electrons and of nuclei that only interact 
amongst themselves. The key assumption in establish- 
ing the MCKS scheme is that local effective potentials 
{Vs,n, Vs,e} exist such that the ground-state densities of 
the auxiliary system reproduce the exact ground-state 



densities {Fq, po} of the fully interacting system. If that 
assumption holds true, the exact ground-state densities 
are given by 



Fo(l) = 

s 



(39) 
(40) 



where x and ipj are solutions of an 7V„-particle nuclear 
and a single-particle electronic Schrodinger equation, re- 
spectively: 

( - E ^ + ^sAW - ) xiM) - (41) 



-—+vsAr)-ee,j ) ^j{r) - 0. (42) 



By virtue of the MCHK theorem applied to the auxil- 
iary system, the effective potentials V5,„(R) and vs,e{^) 
are uniquely determined by the ground-state densities 
{Fq. Po}, once their existence is assumed. They are given 

by ' 



Vs,nm = VF„„(R) 



^-Eu.HxclF^jO] 



(5p(r) 



(5F(R) 



(43) 



(44) 



In this procedure we require the nuclear wavefunction 
X to have the same symmetry-properties under the in- 
terchange of nuclei of the same type as the exact wave- 
function of the interacting system (this will also be re- 
quired along the adiabatic connection to be discussed 
later in the paper). The last terms on the right-hand 
sides of Eqs. (|43fl and (|44|l represent the potentials due 
to all non-trivial interactions of the system, i.e., they 
contain the Hartree-exchange-correlation (Hxc) effects of 
the electron-electron and electron-nuclear interactions as 
well as mass-polarization and Coriolis effects and the in- 
fluence of the true external potentials U. As seen in 
Eqs. H43|) and H44|) . these potentials are given as func- 
tional derivatives of the UHxc energy functional defined 
by 

Su,Hxc[F,p] :- F[F,p] - rs,«[F] - TsAp]- (45) 

This quantity represents the central quantity of the 
MCDFT and contains all many-body effects except the 
purely nuclear correlations. We note that, in the case of 
vanishing external potentials U = 0, the nuclear effective 
potentials Vs_„(R) and the conjugated density, i.e., the 
nuclear density matrix F(R) are invariant under transla- 
tions. Therefore, the nuclear center-of-mass can be sep- 
arated off in Eq. (|41|) . reducing the number of degrees 
of freedom by three. We will illustrate this procedure in 
our applications later. 
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In order to derive the above representations of the ef- 
fective potentials, we consider the energy functional of 
the auxiliary system introduced above: 



+ / dr p(r)ws,e(r)- 



(46) 



As noted before, the nuclear-nuclear interaction Wnn is 
included in the 'external' potential Vs.nC^)- The nonin- 
teracting kinetic-energy functional Ts,e[p\ is the one fa- 
miliar from purely electronic DFT, 



TsAp] 



min($|Te 



(47) 



where the minimization is over all electronic Slater de- 
terminants $ yielding p. Similarly, the nuclear kinetic- 
energy functional is given by 



s,„[r] = min(x|T„|x>- 



(48) 



In contrast to the electronic wavefunction (f>, the nu- 
clear wavefunction x is not a Slater determinant, but 
a correlated many-body wavefunction, since it minimizes 
T„ under the constraint of generating the diagonal of 
the nuclear iV„-particle density matrix. We note that, 
although X is an interacting many-body wavefunction, 
Ts^n is not the interacting nuclear kinetic-energy func- 
tional T„[r,p] = (^'[r,p]|f„|^'[r,p]), since ^T^p] min- 
imizes {^\f + W + U\^) (for given densities {r,p}), 
therefore including all electron-nuclei interactions as well 
as mass-polarization and Coriolis couplings. Assum- 
ing the densities {r,p} to be noninteracting {Vn,Ve}- 
representable, the minimizing states of H47|l and (|48|l . 
i.e., the states minimizing the kinetic energy for given 
{r,p}, are obtained from Eqs. l|IT)l and with the 
potentials uniquely determined by the Euler equations 
following from H46() : 



<5r(E) 

STsAp] 



5p{v) 



+ «s,e(r) 



(49) 
(50) 



Returning to the interacting problem, we decompose the 
functional -^'[r^p] according to Eq. (|45|l . Employing this 
definition in the variational equations H^^7|l and l|38(l of 
the interacting problem and comparing them to the Eu- 
ler equations (|49|l and H50|) . we find that the effective 
potentials which reproduce the exact densities from the 
auxiliary system are indeed given by Eqs. (|43|l and (|44|1 . 

Eqs. (1231- lis, 63) and ^ constitute the MCKS sys- 
tem. Since the effective potentials depend on both den- 
sities, the MCKS equations H41|l and 142|l are coupled, 
reflecting the mutual influence of electrons and nuclei on 
each other, and have to be solved self-consistently. We 



emphasize that Eq. H42() . although similar to the usual 
electronic KS equation, does not parametrically depend 
on the nuclear configuration. Instead, the information 
on the nuclear distribution is already included through 
the functional dependence on T. Considering the nu- 
clear MCKS equation (|41|l . we again realize its simi- 
larity with the conventional nuclear BO equation. Yet, 
no BO approximation has been used to derive Eq. H41|l . 
In contrast, since the MCKS scheme provides the ex- 
act ground state, all non-BO effects are, in principle, 
included. Whether or not the non-BO effects are repro- 
duced in practical applications depends, of course, on the 
quality of the approximations employed for £^u,Hxc[r,p]- 
We also note that in the absence of external potentials the 
potential Kg n,(R) has the same symmetry properties as 
the BO-energy surface under rotations and translations. 



ys,n(OE + a) = Vs.„(R) 



(51) 



The way it will transform under interchange of like nuclei 
will depend on the way we choose the body-flxcd frame. 
It is also important to realize that when solving the nu- 
clear equation (|41|1 we must look for the solution y(Rs) 
with the lowest energy under the constraint that it has 
the correct symmetry under interchange of nuclear space- 
spin coordinates , i.e. the symmetry that was imposed 
by the constrained-search. Like the nuclear BO equation, 
the nuclear equation l|41|) is still a many-body equation. 
Therefore, its solution will, in general, be rather com- 
plicated and further simplifications are highly desirable. 
Typically, one first splits off the nuclear center-of-mass 
motion and the global rotations of the molecule. Then 
the remaining nuclear degrees of freedom are transformed 
to normal coordinates, in terms of which the problem is 
treated in a harmonic approximation, possibly including 
anharmonic effects in a mean-field fashion 36J. However, 
due to the generality of the method, different treatments 
appropriate for different physical situations can be used. 



III. ANALYSIS OF THE FUNCTIONALS 
A. Decomposition of the Energy Functional 

In the last section, the foundations of the MCDFT 
were developed. We derived a formally exact scheme, 
which provides a way to calculate ground-state proper- 
ties of MC systems. For any practical application, the 
functional -Eu.Hxc needs to be approximated. In order 
to gain more insight in the construction of such an ap- 
proximation, this section discusses a number of rigorous 
properties of this functional. 

Following [s^ , we start out by decomposing the UHxc 
energy functional (|45|) in parts associated with its var- 
ious interactions. To this end, we define the following 
quantities: 



F'[p\ := min(^0|re + #ee|^) 



(52) 
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F^"[r,p] := min (^|T„ + re + W^|V') (53) 

7MPc[r, p] ■■= min {xp\fn + + TmpC + W^IV') 

-F-^[r,p] (54) 
C/Hxc[r,p] := F[r,p]- mm {^\f + W\iP). (55) 

The first term represents the electronic functional which, 
by construction, is identical to the functional of 
standard electronic DFT, first introduced in j^j. Usu- 
ally, this quantity is split according to 



ree functioi 
„,P(r)p(r') 



where is the electronic Hartree functional 



E^[P] ■■=\j drdv' 



(56) 



(57) 



and where the electronic exchange-correlation functional 
E'xc is defined by Eq.jSSJ. In contrast to F'^, the second 
functional F^^ also includes the nuclear kinetic energy 
as well as the electron-nuclear interaction, but still ne- 
glects mass-polarization and Coriolis effects and the in- 
fluence of the external potential U. As discussed later 
on, the functional F"^" thus includes in particular the ef- 
fects arising from the electron-nuclear correlation. The 
first term on the right-hand side of Eq. (|54|l addition- 
ally contains the mass-polarization and Coriolis terms 
of the kinetic-energy operator. Therefore, the differ- 
ence between mm{ip\T + W\ip) and F*^" is responsible 
for mass-polarization and Coriolis effects and thus de- 
noted by TmpC- Similarily, the last term denoted by C/hxc 
takes care of all effects introduced by the true external 
potentials U. Consequently, if no true external fields are 
applied to the system, C/hxc vanishes identically. 
Inserting Eqs. into Eq. leads to 



X 



where 



(58) 



(59) 



Eq. H58|) provides a decomposition of the Hxc energy 
functional in its natural contributions. The first part, 
given by E^ and E^^, describes the Coulomb interactions 
among the electrons. It is important to note that these 
functionals are, by construction, identical to the ones 
familiar from standard electronic DFT. The electron- 
electron interaction can therefore be treated in the fa- 
miliar way, namely by using the widely investigated and 
highly successful approximations for the electronic xc en- 
ergy functional E^^[p]. The last term of Eq. ifSSl) was 
constructed to incorporate all effects arising from the 
presence of true external potentials U. As already men- 
tioned above, these terms are not of a single particle form 
in the transformed coordinate system and have to be 



treated similarly to the interaction terms. The functional 
t^Hxc[r,p] provides a means of dealing with these effects. 
Similarly, the fourth term of Eq. I|58|) incorporates all ef- 
fects due to the mass-polarization and Coriolis terms. At 
least for ground-state properties, this term is expected to 
be unimportant and can be neglected in most situations. 
If such effects are, on the other hand, important in a given 
physical situation, they can, in principle, be included in 
the calculation by taking the functional Tmpc[^,p] explic- 
itly into account. Finally, the term E^[r^p] contains all 
effects due to the electron-nuclear interaction. Its analy- 
sis will be continued in the next section. 

The decomposition (|58|l of the energy functional 
£■11,11x0 is obviously not unique. However, the charme 
of the above prescription lies in the fact that, firstly, 
parts like the purely electronic functionals are already 
well known such that one can rely on existing approxi- 
mations for these functional. Secondly, the functionals 
(|53|l - (|55() contain, by their very construction, just the 
effect of one specifically chosen interaction. This, in par- 
ticular, guarantees that the functionals E^^^[p], E^[T^p] 
and T'MPc[r,p] are universal in the sense that they do 
not depend on the external potentials and can therefore 
be employed for all systems independent of the applied 
external fields. All effects arising from the external po- 
tentials are subsumed in [/hxcP.p]. 



B. The electron- nuclear energy functional 

Using the decomposition (|58|l . the well-studied elec- 
tronic Hxc energy functional as well as the - at least 
for ground-state properties - presumably negligible mass- 
polarization and Coriolis contribution were separated off 
in the functional £"11,11x0 • In this section, we discuss the 
functional E^ [^,p] which contains the many-body effects 
due to the electron-nuclear interaction. We will derive an 
equation for this functional that is of a suitable form to 
be used in our approximations later. To do this we will 
use the familiar coupling constant integration technique 
of standard density-functional theory. To begin with, we 
consider the Hamiltonian 



Hx = Tn + Te + Wee + XWen + K,A + ^e,A, 



(60) 



where a non-negative coupling constant A, scaling the 
electron-nuclear interaction, has been introduced. As 
usual, the potential Vx = Vn^x + Ve,x is chosen such that 
the densities remain fixed: Fa = F and px ~ p, indepen- 
dent of the coupling constant A. Employing the coupling- 
constant integration techniq ue 37, 38, .39] adapted to the 
electron- nuclear interaction jsj, the electron- nuclear He 
energy functional H59|) is rewritten as 



mi 
min 



= min (*^|f„+fe + W^ee+ Allien 1*^) 



- min (*^|f„ + fe + W^ee + AT4^e„|^^; 



A=0 
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£ dX ^ (Vl/-;;^^|f:„ + fe + Wee + XWenl^'P^;'^) 



min,A\ 



(61) 



where \]>™'^'^ denotes the minimizing state of (r„ + 3^6 + 

Wee + ^Wen) generating the given densities {T,p}, and 
a Hellmann-Feynman-type theorem was used in the last 
step. Therefore, the electron-nuclear energy functional 
is given by 

EITA^A (62) 
= J d~"Rr(E) J drM/e4E,r)7"""[r,p](r|E), 
where 



|7^-lr-R„ + RcA^Jv| 



(63) 



|r~7^(R„ -Rcj\/A')| 
The electronic conditional density 7 is defined by 

7(r|E) ■■=NeY. f d'^'-'r |*(i^, Ei)l Vr(E), (64) 

and 7 represents the coupling-constant average of 7. The 
conditional density satisfies important sumrules that we 
will use later for the construction of approximate func- 
tionals: 



Ne 



P(r) 



j dv 7(r|R) Va (65) 



dRr(R)7(r|R) (66) 



By virtue of Eq. , the He energy can be interpreted 
as the electrostatic interaction energy of the (coupling- 
constant averaged) electronic density for a fixed nuclear 
configuration with the point charges of the correspond- 
ing nuclei, averaged over the nuclear distribution. This 
interpretation will play an important role in our later 
development of approximate functionals. 

In order to gain further insight into the electron- 
nuclear He energy functional, we now establish a connec- 
tion between the MCDFT scheme and the conventional 
BO method which provides a highly successful treatment 
of electron-nuclear correlation. To that end, we decom- 
pose the total wavefunction into an adiabatic product 
according to 



*r;p"(m,Rs) = x(R^) S[r,p](m|Rs), 



(67) 



where x is the nuclear wavefunction generating the nu- 
clear density matrix T and S[r_p] is an electronic state 
normalized to one for every nuclear configuration Rs : 



j d^=r|s[rp](ig|Ei)P-i. 



(68) 



We note that the decomposition 1)6 7|) is actually an ex- 
act representation of the correlated electron-nuclear wave 
function ^ El SIS El El and that the factors x 
and S are unique l4g| up to within an Rs-dependent 
phasefactor. However, it is important to note that the 
electronic state S is not identical to the usual electronic 
BO state. Even if non-BO effects were neglected, S 
would not be identical to the electronic BO state since 
X and S are required to reproduce a given set of den- 
sities (r,p) (the two electronic wavefunctions only be- 
come equivalent, if S[r,p] is evaluated at the BO den- 

according to 



sities (r,p)=(r^o,pB°)). Instead, S[r^p] is expanded 



(69) 



where } denotes a complete set of (BO) eigenfunc- 

tions corresponding to the electronic (clamped-nuclei) 
Hamiltonian He := T; + Wee + Wen- Employing Eq. 
together with H68|l and assuming that S[r, p] is real, 
the electron-nuclear He energy functional is given from 
Eq.(|nni) by 

i?Hc[r,p] (70) 

= Y,[ d^"R |x(Ei)P(s[r,p]|f„ + i?e|s[r,p])e - f-^M 

where the index "e" at the bracket indicates that the inte- 
gration is over electronic coordinates only. Using Eq. I|69(l 
we then obtain 

EZ[^A (71) 

= Y,[ d^"R \xm\^ { l«fe[r,p](Ei)i' • er(E)'5, 

s •' k,l 

+ 4[r.p](Ei) (snT„|sfO)ea4r,p](E|)} - F^[p\ 

where 



k,l 



6r(R) = (sriifeiHr) 



"BO I 



|"BO\ 



(72) 



represents the fcth BO potential-energy surface (PES) 
and the index "e" at the bracket indicates that the inte- 
gration is over electronic coordinates only. On the basis 
of Eq. H71|l . one can interpret as the potential en- 
ergy of the nuclei, where the nuclear distribution lies in 
a potential hypersurface, which is composed of adiabatic 
BO-PES weighted with the coefficients a^p^p] as well as 
nonadiabatic corrections to it. Of course, the coefficients 
Qk and their functional dependence on the set of densities 
[r^p] is unknown at this point. However, Eq. (|7T|l helps us 
in gaining a better understanding of the electron-nuclear 
He energy functional and establishes an - at least - for- 
mal link to the BO scheme which is further exploited 
when the effective potentials are discussed later on. 
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C. Concerning the true external potentials 

Similar to the techniques employed in the last section, 
the coupling-constant integration can be employed to de- 
rive an expression for the functional ?7Hxc[r,p] which sub- 
sumes the many-body effects arising from the true exter- 
nal field. In analogy to Eq. (|62|l . one obtains 



(73) 



= J d^"R r(E) J dr C/e,t(E, r) 7™ p](r|E) , 

where 7""" again denotes the coupling-constant average 
with respect to coupling constant fi of the conditional 
density 7"""'^''[r,p] corresponding to the states vp™'" 

minimizing {T + W + fiU) . We further defined [/cxt (R, r) 
to be 



f/cxt(E,r) := ^C/oxt.n(E) 



-Ucxt. 



,in-'r+RcMN) (74) 



It has to be noted that the conditional densities appear- 
ing in Eqs. H62|l and (|73|1 are not identical since the cor- 
responding states vji™'"^^ and ^'p"p minimize different 
expressions. This is a direct consequence of the defini- 
tions chosen in Eqs. 1)52(1 - (|55(l . In particular, this choice 
guarantees that E^^^, Tmpc, and E^l are independent of 
the true external potential Ucxt, i-e., these functionals are 
universal; all effects stemming from Ucxt are contained 
exclusively in the functional J7hxc- 

By virtue of the above discussion, the influence of the 
true external potential has to be treated similar to an 
interaction. As already mentioned above, this compli- 
cation is an immediate consequence of the necessity to 
transform to an internal reference system for the formu- 
lation of the MCDFT scheme. Of course, in the numer- 
ous cases discussing the properties of isolated systems, U 
vanishes and the MCDFT formalism reduces to the one 
given in 's^. If, on the other hand, a true external po- 
tential is applied to the system, approximations for ?7hxc 
are needed. In the simplest case, the electronic condi- 
tional density 7 is replaced by the electronic density p, 
leading to a Hartree-type approximation for C/hxc- Such 
an approximation will be especially valid in the case of 
well localized nuclei, as discussed later on. 



D. Analysis of the Effective Potentials 

In the last sections, the Hxc energy functional of 
Ea. ((^ was discussed. According to Eqs. (gSJ and ((I^ . 
this quantity gives rise to the many-body contributions 
of the effective MCKS potentials. Explicitly, the UHxc 
potentials are given by 



Vu.Hxc[r,p](E) = 

i'u,Hxc[r,p](r) = 



3Eij,-tixc[^^p] 
6p{r) 



(75) 
(76) 



Employing Eq. (|58|l . the potentials can also be decom- 
posed into the parts associated with the different inter- 
actions, yielding 



^u,Hxc[r,p](E) 

wu,Hxc[r,p](r) 



*'Hxc 



r,pm + ^S^c[^M) 
i/MPc[r,p](E) (77) 

^;g,,[^p](r) + ^;|^[p](r) + <Jp](r) 
vlll[T,p]iv)+VMPc[T,p]iT^), (78) 



where the various potential terms on the right-hand sides 
of the above equations are defined in analogy to H75() and 
(|76|l : The first terms on the right-hand side of Eqs. H75|) 
and ((76|l represent the influence of the true external po- 
tential U and correspond to the derivatives of J/hxc in 
Eq. (|58() . Since the electron-electron interaction is treated 
employing the well-known Hxc energy functional E^^^ [p] 
from standard electronic DFT, the corresponding poten- 
tials and are also identical to the familiar elec- 
tronic Hartree and xc potentials. Furthermore, as for the 
energy functional, the potentials arising from the mass- 
polarization and Coriolis effects are not expected to con- 
tribute significantly - at least for ground-state properties. 
In the following we will concentrate on the Hxc potentials 
arising from the electron- nuclear energy functional E^. 

To start with, we consider the nuclear He potential, 
defined by 



^Hc[r,p](E) 



SEIPA^,p] 

(5r(E) ■ 



(79) 



Employing the representation of i?^" in terms of the 
coupling-constant averaged conditional density, Eq. l(B^ . 
the nuclear potential can be split in two parts. 



vi-[T,pm) = v^-j^MR) + KTsp[r,p](E) (80) 



where 



K?nd[r,p](E) J dr Wen{R,r) ^'"[r, p](r|E) (81) 

is the electrostatic potential due to the electronic condi- 
tional density and 

KTspir.pKE) := J d^"R'r(E') 

xJdrWUR,r) ^r(E) ~ 



(82) 



defines a response-type contribution to the electron- 
nuclear correlation potential. We note that the condi- 
tional potential Vj,™ completely determines the electron- 
nuclear He energy: 

i?^"jr,p] = J d^"Rr(E)v;ond[r,p](E)- (83) 

In the following analysis we restrict ourselves to situa- 
tions where the full electron-nuclear wave function 'I'p"" 
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can be factorized into a nuclear spin function times a re- 
mainder not depending on s. This is exactly true, e.g., 
for diatomic molecules or when the nuclei are spin-zero 
bosons. In many other cases, this factorization represents 
a good approximation. Under these circumstances, the 
wave function S in Eq. (|67(l can be chosen to be indepen- 
dent of s, and likewise the expansion coefficients in 
Eq.lEHirso that Eq.lO reduces to 

i?Hc[r,p] (84) 
= / d^"R r(E) { WiTMRf ■ ^nmsk,i 

k,l 

+ 4[rp](E)(Snf„|SfO),aarp](^)}-F^[p], is given by: 

Comparing Eq.ljSSJl and (|84|l . the conditional potential 
(|81|1 can be expressed in terms of the BO-PES: 

KZA^Mm = E l«^-[r.p](E)r^r®) 

k 

+ E4[r,p](E) (s^°|f„|sr)ea4r,p](E) 



k,l 



(85) 



This equation provides a useful tool to interpret the effec- 
tive nuclear MCKS potential: The first term in H85() is a 
weighted sum over different adiabatic BO-PES, whereas 
the second one describes adiabatic and nonadiabatic cor- 
rections to it. The last term in equation H85|l . F'^[p\, just 
yields a constant shift and is included in the potential 
to maintain the same zero-energy level within the BO 
and MCKS schemes. Considering the case that the BO 
approximation accurately describes a specific system, we 
realize that, in the first sum, only the lowest coefficient 
ao survives and the second sum is negligible, provided 
the potential is evaluated at the ground-state densities. 
Therefore, Fif^(E) ~ eo°(E), and the nuclear MCKS 
equation reduces to the nuclear BO equation in the limit 
considered here. We emphasize, however, that the way 
to evaluate this potential differs in the MCKS and BO 
methods. Whereas, in the latter, an electronic equa- 
tion has to be solved for each nuclear configuration, the 
MCKS potential is determined by the functional deriva- 
tive 5E^^\r^p\/5T. Inserting the ground-state densities 
then yields a potential which, as a function of R, is very 
close to the BO potential (in the case discussed here). 
Furthermore, we can conclude that the response part of 
the nuclear potential, Eq. (|S^ . has negligible influence 
for such systems. If, on the other hand, nonadiabatic 
effects - e.g. close to level crossings ~ are encountered, 
the coefficients Ok in Eq. H85() will, as a function of the 
nuclear configuration, achieve a natural diabatization. 
One should also note that the electronic wavefunction 
S is, in general, complex at points of degeneracy. There- 
fore, one obtains another contribution to the nuclear po- 
tential, which is responsible for Berry-phase effects |47l | 
|60j| . In addition, the response part of the nuclear poten- 
tial might contribute appreciably. In summary, Eq. H85|) 



shows that the (exact) nuclear effective potential reduces 
to the lowest-energy BO-PES, if nonadiabatic contribu- 
tions can be neglected, but also contains in principle all 
non-BO effects. Whether or not they can be recovered 
in an actual application crucially depends, of course, on 
the level of sophistication of the approximation used for 

Employing again Eq. 1)62(1 , the electronic potential due 
to the electron-nuclear interaction, defined by 



^Hc[r,p](r) := 



Sp{r) 



(86) 



y^"jr,p](r) = J d^"Rr(E) 



^'w ,/7'"'"[r,p](r'|R) 

X / dr'Wen{R,r') — - — (87) 

— dp{r) 

This expression appears rather complicated to evalu- 
ate. If, however, the nuclear probability distribution 
is sharply peaked around an equilibrium geometry R , 
only configurations around R^^ will substantially con- 
tribute to the above integral. Then, the calculation of 
the electronic Hxc potential simplifies to 



<[r.p](r) 



d^^"Rr(R)VFe„(R,r). 



(88) 



This potential represents the electrostatic (Hartree) po- 
tential due to the nuclear charge distribution acting on 
the electrons. Since the nuclear ground-state densities of 
many molecules are indeed strongly localized functions - 
in other words: the nuclei behave almost classically - we 
expect the Hartree approximation for the electronic po- 
tential to be sufficiently accurate for such systems. If, on 
the other hand, the assumption of nicely localized nuclear 
densities breaks down, one needs to incorporate correla- 
tion contribution to vfH arising from the electron-nuclear 
interaction. 



E. The Limit of Classical (Point-Like) Nuclei 

In this section, we investigate the limit of classical, i.e., 
perfectly localized nuclei. Assuming identical zero spin 
nuclei for the ease of notation, the nuclear density matrix 
reads 



r'='^^^(Ri...R. 



'JV„ 



^ E n (R-P(a) - R-a.o) 



(89) 



V a 



where the sum is over all iV„! permutations of the nu- 
clear coordinates and R^ denotes the positions where the 
nuclei are located. Note that by this classical form of 
the density matrix we have broken the translational and 
rotational symmetry of the density matrix as presented 
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in Eg. pifl . In the following we investigate the conse- 
quences are of this form for the diagonal of the nuclear 
density matrix. First, we consider the electronic density. 
In terms of the coupling-constant dependent conditional 
density, it is given by 



p{r)= I d^"Rr(R)7'^(r|R) 



(90) 



We recall that p{r) does not depend on the coupling con- 
stant A, since the external potentials are chosen such that 
the densities remain unchanged. Inserting Eq. H89|> into 
(EOI yields 



(91) 



i.e., in the limit of classical nuclei, the electronic density is 
identical to the conditional density evaluated at the posi- 
tions of the classical nuclei. This quantity, in fact, serves 
as the basic variable of standard electronic DFT employ- 
ing the BO approximation: p°^"'"'^'-'(r) = 7(r|R,p). We 
therefore conclude that the MCDFT presented here re- 
duces to the standard formulation of DFT in the limit of 
classical nuclei. 

Inserting Eq. Ip^ into Eqs. and JTSJ), we readily 
obtain the expressions for E^l and J/hxc in the classical 
limit: 



;lass 



,p] = / drp(r)M^e„(Eo,r) 



dr p(r)[/ext(Eo,r) 



(92) 
(93) 



Thus, in the limit of classical nuclei, the Hxc energy func- 
tional reduce to the classical electrostatic (Hartree) in- 
teractions and correlation contributions vanish [s^ ]. 

The corresponding electronic potentials, following 
from Eqs. (|H2J) and then read 



en [T^class 
^Hc 



•'Hxc 



We 



C/cxt(Eo,r). 



(94) 
(95) 



The first quantity is identical to the classical Coulomb 
field of the nuclei, whereas the second one describes the 
influence of potentials applied externally to the system. 
Both quantites together represent the "external poten- 
tial" in BO-based DFT, reflecting again its coincidence 
with the MCDFT in the limit of classical nuclei. From 
this perspective, one also might consider Eq. (|88|l as the 
natural extension of Eq. (|94|l to nuclear distributions 
which are localized but still exhibit a flnite width. 



IV. APPLICATIONS 
A. Diatomic molecules 

Having discussed the foundation and some formal 
properties of the MCDFT, we proceed with the appli- 
cation of the theory to the case of isolated diatomic 
molecules. 



However, as mentioned above, we treatment of the Nn- 
body nuclear MCKS equation has to be discussed prior to 
actual applications: Since "true" external potentials are 
absent, i.e., U = the system is translationally invariant. 
Accordingly, the nuclear MCKS potential is required to 
behave as Vs.n{'R-i,'R-2) = V5,„(Ri — R2). The nuclear 
equation has then the form 



2Mi ^1 



1 



2M2 



R2 



+Fs,„(Ri - R2) - e„jx(RiSi,R2S2) = (96) 

Then the nuclear CM motion can be separated off and the 
problem can be reformulated in terms of the internuclear 
separation R := R2 — Ri. The nuclear function x has 
the general form 

x(RiSi,R2S2) =?7(Rcm)C(Ri -R2)e(si,S2) (97) 

where 77(Rcj\/) is a plane wave state depending on 
the center-of-mass coordinate Rcm — (AfiRi + 
A/2R2)/(Mi + M2) and is expHcitly given by 



?7(Rcm) = 

vV 



k-Rr 



(98) 



where V is the total volume of the system. The rela- 
tive wavefunction ^(Ri — R2) satisfies the nuclear MCKS 
equation 



VsA^) - e„ C(R) = 0, 



(99) 



where pn = M1M2/ {Mi + M2) denotes the reduced nu- 
clear mass. The function is a nuclear spin function 
depending on the type of nuclear species. For instance, 
for the H2 with two protons as nuclei the nuclear ground 
state is that of para-hydrogen where the function 9 is an 
antisymmetric spin function and consequently the func- 
tion must be even, i.e. ^(R) = ^(— R), to preserve 
overall anti-symmetry of the wavefunction under the in- 
terchange of the two protons. For the density matrix F, 
one obtains after integrating out the spin function 

r(Ri,R2) = -i|e(Ri - R2)P = -^r(R). (100) 

Therefore, the diagonal of the nuclear density matrix, 
which we often - and somewhat unprecisely - refer to 
as "nuclear density" is indeed a single-particle quantity 
describing the probability of finding the two nuclei sepa- 
rated by R. 

It remains to discuss the electronic coordinates. For our 
diatomic molecule we determine the Euler angles by the 
requirement that the internuclear axis be parallel to the 
z-axis in the body-fixed frame, i.e. 7?.(Ri — R2) = Rbz 
, where R = |Ri — R2I. For the special case of the 
diatomic molecule only two Euler angles are needed to 
specify the rotation matrix TZ. The electronic coordinates 
in the body- fixed frame are then obtained using Eq.®. 
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With this transformation the electron-nuclear attraction 
attains the form 



WenCR, r) 



Zi 



M2 



■Re, 



Ml 

Mnuc 



Re, 



(101) 



As for any other DFT, explicit approximations have to 
be employed for the energy functional i?u,Hxc[r,p] of 
Ea. (|58|l . Since no "true" external potentials are present 
in the case discussed here, C/Hxc[r,/o] vanishes identically. 
Furthermore, following Sec. IIII Al the purely electronic 
part of the energy functional can be treated by using 
the familiar approximations for the electronic xc energy 
functional E^^[p\. For all systems containing more than 
one electron, we will employ the well-known LDA ap- 
proximation. We emphasize that it is not the purpose 
of this work to investigate new approximations for the 
electronic xc energy functional. Instead, we aim at an 
analysis of the previously not much studied He energy 
functional arising from the electron-nuclear interaction. 
To that end, we restrict ourselves to work within the 
LDA approximation for E^^.[p]; different - and more so- 
phisticated - approximations for the electronic xc energy 
functional would only result in some minor quantitative 
changes of the analysis presented below. Furthermore, 
we note that the mass-polarization and Coriolis effects 
are not expected to contribute substantially to ground- 
state properties. Therefore, only the diagonal part of the 
mass-polarization term which leads to a reduced mass 
in the electronic MCKS equation is accounted for and 
the remaining parts are neglected in all practical calcula- 
tions. With these assumptions the electronic Kohn-Sham 
equations for our problem then attain the form 



where 



-—+vsA^,P]i^) - ^e.j ) = (102) 



numerical integration of the full many-body Schrodinger 
equation. In this way the exact nuclear potential Vs^n for 
this problem was obtained For this case we found 

the exact Vs.n to be almost identical to the BO poten- 
tial except for the case when the nuclear masses were 
taken to be artificially small. This illustrates the point 
discussed before: in situations where the BO approxima- 
tion works well, the nuclear potential Vs^n will be close 
to the BO potential. For this reason the BO potentials 
of the H2 and molecules that we will study below 
will be a good reference to test our approximations for 
the electron-nuclear correlation functional. Of course, 
having obtained a good approximate functional its main 
field of applicability will be cases where the BO approx- 
imation does not work well. 

To conclude these introductory remarks, the numerical 
implementation of the MCKS equations H99() and H102() 
is briefly described. Since Coriolis effects are neglected, 
the effective nuclear potential is spherically symmetric 
and the angular part can be treated analytically. The 
remaining radial nuclear MCKS equation is numerically 
solved on a one-dimensional grid. Furthermore, we ob- 
serve that the z-component of the electronic angular mo- 
mentum is a conserved quantity. Hence, the electronic 
MCKS equation is rewritten in terms of cylindrical co- 
ordinates. For axial-symmetric electronic potentials, the 
angular part can be integrated out and we are left with 
a two-dimensional problem. The resulting Hamiltonian 
is then discretized on a uniform rectangular grid and nu- 
merically diagonalized by employing the Lanczos algo- 
rithm |43|. Due to the use of finite uniform grids, the 
regions around the nuclei are not sampled with high ac- 
curacy, leading to typical discretization errors of about 
0.1% for the systems discussed later in this section. Both 
the nuclear and the electronic equation are solved simul- 
taneously until self-consistency is achieved. 



ws,e[r,p](r) = 



«Sc[r,p](r) 



-t;^[p](r) + <Jp](r) (103) 



and /ie = (Mi + Al2)/{Mi + M2 + l)- The approximations 
used for will be discussed in the subsequent sections. 
Within the same assumptions the effective potential in 
the nuclear equations will be of the form 

VsA^, P](R) = WeeiR) + Vi^[T, p](R) (104) 

where I^ee(R) = Z1Z2/R. 

We finally note that the diatomic molecule is a partic- 
ularly convenient case for studying the nuclear effective 
potential Vs.n since for this case for a given density r(R) 
the potential Vs.n is easily constructed from inversion of 



Fs,„(R) - e„ 



1 V2,yf(R) 



(105) 



2m« v'r(R) 

dimensional mc 
for which the exact nuclear density F can be obtained by 



We have done this for a one-dimensional model of the 



B. The Hartree Approximation for the 
electron-nuclear energy functional 

It remains to find explicit approximations for the 
electron-nuclear He energy functional Ef[l:[T.p]. In the 
simplest case, the electron-nuclear interaction is approx- 
imated by the Hartree energy functional, defined by 

^&"[r,p] J dRidR2drTiR)WeniR,v)p{r) 

= J dRdrTiR)W,niR,r)p{r). (106) 

where in the second step we changed to relative R and 
center-of-mass Rcm coordinates and performed the inte- 
gration over RcAi which eliminated the inverse volume 
prefactor in Ea. l|100|l . By virtue of Eq. (|106|l . the He 
energy functional E^^[r p] is thus replaced by the classi- 
cal electrostatic interactions of the corresponding charge 
distributions and correlation contributions are neglected. 
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BO 


Hartree OAO 


SAO 


-Eo 


1.130 


1.121 


1.122 


1.124 


Ts 




1.069 


1.063 


1.049 


^Hxc,e 




0.627 


0.625 


0.623 






3.496 


3.487 


3.471 


Wnn 




0.679 


0.676 


0.676 


{R) 


1.49 


1.48 


1.49 


1.50 




2.25 


2.21 


2.24 


2.28 


a;[cm~^l 


4137 


7945 


7047 


4282 



TABLE I: Summary of results for the H2 molecule obtained 
from self-consistent solutions of the MCKS scheme employing 
various approximations. For comparison, results from BO 
calculations are added. The electronic interaction is treated 
Within the xcLDA. All numbers (except u>) in atomic units. 





BO 


Hartree OAO 


SAO 


— En 


0.598 


0.591 


0.595 


0.581 


Ts 




0.591 


0.583 


0.574 


rpen 




1.673 


1.662 


1.642 


Wn„ 




0.491 


0.485 


0.487 


(R) 


2.07 


2.05 


2.08 


2.08 




4.30 


4.22 


4.37 


4.39 




2297 


5191 


3248 


2232 



TABLE II: Summary of results for the H2 molecule obtained 
from self- consistent solutions of the MCKS scheme employing 
various approximations. For comparison, results from BO cal- 
culations are added. All numbers (except lo) in atomic units. 



Evidently Eq. (|106|l can also be derived from a product 
(mean-field) ansatz for the electron-nuclear part of the 
total wavefunction. In fact, such a mean- field description 
of the electron-nuclear interaction has been proposed in 
[i^lsoj to study the protonic structure of molecules. 

From the Hartree-energy functional (|1U6|) . the corre- 
sponding potentials, defined in ifTi^ and l^^. are readily 
calculated: 

= J dvWen{R,r)p{r) (107) 

v^(r) = y"dRV7e„(R,r)r(R). (108) 

We note that, within the Hartree approximation, the 
nuclear response potential (|82|l vanishes and the con- 
ditional potential H81|l is given by H107|l . With these 
Hartree potentials inserted in expressions H103() and 
^U^, the MCKS equations (EHJ and ifTU^ are solved 
self-consistently as described above. 

1 . Results 
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FIG. 1: Effective nuclear potential Vs,n{R) for the 
H2 molecule obtained from self- consistent solutions of the 
MCKS scheme employing various approximations. For com- 
parison, results from exact and BO calculations are added. In 
atomic units. 



In the following, the application of the Hartree approx- 
imation to the H2 and H2 molecules is discussed. 

In tables and |n] a selection of results is presented. 
Since the BO approach provides an excellent approxima- 
tion for the system under consideration, we also added, 
for comparison, the results obtained from the BO cal- 
culation (employing the very same numerical procedure 
discribed above). The ground-state results are found to 
be surprisingly good for both molecules. Compared to 
the BO results listed in the first column of the tables, the 
total ground-state energy and the geometry, represented 
by the mean internuclear distance (i?), are reproduced 
up to within an accuracy of about 1% by the Hartree 
method. However, turning towards the harmonic con- 
stants, we realize that the Hartree result is off by about 



a factor of two for both H2 and H2 . Since uj measures 
the curvature of the effective nuclear potentials at the 
minimum, we may expect larger deviations in this quan- 
tity. In fact, this is confirmed by Fig.Oand Fig.^ where 
the effective nuclear MCKS potential is plotted for the 
H2 and H2 molecule. Clearly, the Hartree potential is 
satisfactorily only in a small region around the minimum 
of the potential - which, however, is sufficient for good 
results for the total ground-state energy or geometry. For 
larger R, the potential grows much to fast such that the 
depth of the potentials is largly overestimated. 

In fact from the large R behavior of Ea. (|107|) we see 
that the asymptotic behavior of the effective potential in 
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FIG. 2: Radial nuclear density 47r_R^r(_R) for the H2 molecule 
obtained from self-consistent solutions of the MCKS scheme 
employing various approximations. For comparison, results 
from exact and BO calculations are added. In atomic units. 
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FIG. 3: Effective nuclear potential Vs,n{R) for the H2 
molecule obtained from self-consistent solutions of the MCKS 
scheme employing various approximations. For comparison, 
results from a BO calculations is added. In atomic units. 

the Hartree approximation is given by 

VsA^) {Z1Z2 A^eAf„uc(^ + ^))\ (109) 

whereas the BO potential approaches a finite value in 
this limit. As a consequence of the steep rise of the 
Hartree potential, the corresponding nuclear densities of 
both molecules, shown in Fig|21and Fig 0] are much more 
localized than the BO ones and reflect the wrong shape 
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FIG. 4: Radial nuclear density 47ri?^r(7?) for the II2 molecule 
obtained from self-consistent solutions of the MCKS scheme 
employing various approximations. For comparison, results 
from exact and BO calculations are added. In atomic units. 



of the effective nuclear potential. 

To understand the origin of these deviations in more 
detail, we reconsider the expression for the Hxc energy 
functional, Eq. (|62|l . A comparison of this equation 
with the Hartree energy functional Ij 106(1 shows that the 
Hartree approximation for the electronic conditional den- 
sity is 

7H(r|R) = p(r) VR, (110) 

i.e., the conditional density is independent of R in the 
Hartree approach. In the case that the nuclear density 
is a well-localized function, the approximation H110|l is 
justified for R « Req, (where Rg^ denotes the equilib- 
rium separation) since many quantities of interest then 
only depend on intcrnuclcar separations close to Reg. 
Thus the approximation (|110|l leads to reasonable re- 
sults, as reported above. However, Eq. Hll()|l fails for 
large |R — Req|. This is illustrated in Fig. O where we 
sketched a typical behavior of the electronic density p(r) 
and the conditional density 7(r|R) for |R| « 6 a.u. As 
a consequence, the Hartree approximation cannot be ex- 
pected to be accurate for R > |Reg|, explaining the de- 
viations discussed in Sec. IIVBI 

To summarize, the Hartree approximation provides 
a fair estimate for ground-state properties of diatomic 
molecules, such as the total ground-state energy or the 
equilibrium geometry. However, the nuclear Hartree po- 
tential only reproduces the position of the minimum but 
fails to correctly describe the shape of the exact MCKS 
potential. If one is interested in quantities depending 
more sensibly on the shape of the nuclear potential, the 
Hartree method thus needs to be improved. 
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FIG. 5: Typical behavior of the electronic density p{z) and the 
conditional density j{z\R) ~ 6 a.u.) for a diatomic molecule, 
plotted along the internuclear (z) axis. In atomic units. 

C. An Approach Based on Atomic Orbitals 

From the analysis of the preceding sections, we were 
lead to the conclusion that one needs to go beyond the 
Hartree approximation for the electron-nuclear energy 
functional E^J^T^p]. Moreover, in view of the fact that 
the deviations between the Hartree and the exact (BO) 
potentials are rather large, the Hartree potential is not 
necessarily a good starting point for approximations - in 
contrast to the situation for the electron-electron inter- 
action. This is particularly true for properties which do 
not only depend on internuclear distances close to Req. 
For instance, to go beyond the Hartree approximation, it 
might be tempting (as usually done for electronic corre- 
lations in standard DFT) to define a "hole density" 

7c(r|R) := 7(r|R) - p(r) (111) 

which measures the deviations of the density from the 
conditional density. Inserting this definition into Eq. I|62|) 
then leads to partitioning of the Hartree-correlation func- 
tional into a Hartree and a correlation part. However, as 
displayed in Fig.© the hole function 7c(r|R) is large 
for almost all R (except R w Roq) and has to account 
for large correlation corrections. Therefore at the present 
point, it appears most promising to approximate the con- 
ditional density directly. To that end, we first recall that 
7(r|R) represents the probability of finding an electron at 
r, provided the nuclei are separated by R. If we would 
go back to a BO description such a quantity would be 
naturally described by the electronic density calculated 
within the BO approach as described in Eq. (|91|) in which 
case we have 

7^°(r|R) = f;|Ci:^(r)p (112) 

where ^^'^ are the electronic single-particle orbitals for 
internuclear distance R as in standard density functional 



theory. In order to obtain Ea. H112() in the BO-limit we 
approximate the conditional density 7 in the spirit of a 
linear combination of atomic orbitals (LCAO) approach, 
i.e., 

7(r|R)«5]7,(r|R) (113) 

i=i 

where 

7,(r|R) - \<t^f{vA) + <I>^{vb)\^ , (114) 

where (j3^^^{r) denotes an atomic-type orbital and the 
factor i'j(R) is included to ensure normalization of the 
orbital conditional density, i.e. 

1 = J dr7, (r|R). (115) 

This ensures normalization of the electronic Kohn-Sham 
orbitals in the BO limit. Instead of following the stan- 
dard LCAO approach where 7(r|R) is constructed from 

given atomic orbitals {0j }, we will determine suitably 
"optimized" atomic orbitals (OAO) from given densities 
r and p, i.e., the atomic orbitals are represented as func- 
tionals oi the densities: 0^^^(r) = [T, p]{r). Insert- 
ing these orbitals in (|114|l then leads to an approximation 
of 7[r, p] and therefore, by virtue of Eq. (|62|l . to an ap- 
proximation of E^l[r, p] as functionals of the densities F 
and p. 

In order to find the OAO, we first note that, given 
atomic orbitals 4>f{Y) and 0^(r), normalized bonding 
and antibonding molecular orbitals for a fixed internu- 
clear distance R can be obtained from 

where SjiR) :— J dr (j)^ {r A)4'f {y b) denotes the overlap 
integral and the atomic orbitals are assumed to be real. 
For the purpose of the section, we are concerned with the 
reverse problem: Given molecular MCKS orbitals, how 
can we construct corresponding atomic orbitals? 

As the crucial idea, we identify the electronic MCKS 
orbitals {iy9j(r)}, i.e., the solutions of the electronic 
MCKS equation (|102|l . with the bonding and antibonding 
orbitals of Eq. H116|l evaluated at the mean internuclear 
distance {R): 

^,(r) = ^f^)^,(r) (117) 

where Lp-j denotes the antibonding counterpart of ipj . Us- 
ing Eqs. H116|) and H117|) . we can solve for the atomic 
orbitals, yielding 
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</.f[r,{^,}](r) = 
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M2 




Mnuc 






Ml 


"Pi 1 




-^'-^nuc 






Ml 


V'J 1 


(- 





{R)e, 119) 



(i?)e,jl.20) 
(121) 

The value of the overlap Sj = Sj{{R)) is not deter- 
mined by the above procedure and has to be supplied 
additionally. A simple estimate is obtained from the 
overlap of unperturbed (hydrogenic) orbitals. Using this 
prescription, the atomic orbitals can be calculated from 
Eqs. H119() and (I120() . We therefore determined the 
atomic orbitals as functionals of the nuclear density and 
the electronic MCKS orbitals or, by virtue of the MCHK 
theorem, as implicit functionals of the densities F and p. 

For the case of homonuclear diatomic molecules, the 
above equations yield (f>f'{r) — (f)^ {~r). We note that the 
atomic orbitals are not required to have a definite symme- 
try with respect to parity transformations. As a matter 
of fact, we do not expect them to be symmetric; instead, 
one may view them as orbitals which are centered on one 
nucleus and polarized by the presence of the second nu- 
cleus. Of course, for homonuclear molecules, the linear 
combinations (|116|l are properly symmetrized, i.e., the 
molecular orbitals can be classified either as gerade or as 
ungcrade states. 

Employing Eqs. H119|l and (|120|l . we readily obtain an 
approximation of the conditional density as an (implicit) 
functional of the densities is obtained: 



70^°[F,{^,}](r|R) 



(122) 



TVe 

E 



/>/[F,{^,}](r^) 



0f[F,{^,}](rs) 



In the asymptotic R ^ 00 regime, Eq. (|122|) reduces to 
the correct asymptotic form, i.e., the conditional den- 
sity is given by the sum of two atomic densities How- 
ever, for the self-consistent ground-state solution of the 
MCKS scheme employing Eq. 1)12211 . the atomic orbitals 
are not the unperturbed orbitals representing the ground 
state of the dissociated fragments, but rather polarized 
orbitals which are optimized for the molecular ground 
state. Therefore, even for i? — > 00, the conditional den- 
sity j'^^'^ is not exact, although it should improve on the 
Hartree behavior. 

Employing Eq. (|122|l . we obtain an expression for the 
electron-nuclear interaction energy. 



w: 



OAO 

en 



r,p]= J dRF(R 
J drWen{R,r) 7°^' 



(123) 



F,{^,[p]}](r|R). 



In order to calculate the He energy functional from this 
expression we could perform the coupling-constant inte- 
gration by using an approximation for the A-dependence 



of the conditional density in Ea. l61|) . Alternatively, we 
employ Eq. ifS^ . yielding 



TTiOAO 



[r,p] 



dR F(R) 



(124) 



J dr I^e„(R,r) 70A0(r|R) + F^[70A0](r) _ p^^p^^ 

where nonadiabatic terms have been neglected and F"^ 
is the universal electronic functional defined in H52|l . 
Eq. (|124|l represents the central result of this section. 
In order to use this approximation in a self-consistent 
MCKS calculation, the effective potentials have to be 
calculated. The nuclear conditional potential is readily 
evaluated, yielding 



yOAO 
•^cond 



(R) 



T/OAO (Ty-\ 
'^cond,Wl-'^7 



_ T/OAO 

' r 



cond,T(R) 



OAO 
condjXC 



where 



v: 



GAG 
cond,W 



(R) = J drW^e„(R,r)7°^°(r|R) 



(125) 



(126) 



taOAO 
''cond.T 



(R) 



1 



f |V7}^^"(r|R)p 

"-^ ^,OAOc 



TsAp] 

1 



dvdv' ^ 



7^-^(r|R) 

(127) 

o^O (r|R)70^0(r^|R) 



|r-r'| 



-rrOAG 
''^cond.xc V 



2 

^hW (128) 

,(R) = i?^j7°^°](R)-i?^c[p]- (129) 

The first term on the right-hand side of Eq. p25|l rep- 
resents the functional derivative of VFg^'^^p, with re- 
spect to the nuclear density. The remaining terms are 
responsible for the i?-dependence of the electronic con- 
tributions to the Hxc energy functional. The functional 
derivative of Eq. (|124|l is defined up to an arbitrary con- 
stant. For this reason the terms H127|l - H129l) are defined 
such that they vanish when 7(r|R) w p(r). The nuclear 
response potential is neglected for the reasons explained 
already above. In order to calculate the electronic poten- 
tial from (|124|l . one would have to resort to the optimized 
effective potential (OEP) method [sil Is^. 53], since the 
OAO energy functional depends explicitly on the elec- 
tronic MCKS orbitals {fj} and therefore implicitly on 
the electronic density. This, however, would lead to a 
rather complicated integral equation. On the other hand, 
it was shown in Sec. lIVBl that the electronic Hartree po- 
tential is sufficiently accurate for the systems considered 
here. Therefore, v^''{y) will be used as an approximation 
for the electronic He potential in the current context, 
too. Having derived the effective potentials, the MCKS 
scheme is solved self-consistently. 



1 . Results 

In the following, the results of the OAO appoach are 
presented for the H2 and the H2 molecules. First of all. 
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FIG. 6: Atomic orbital obtained for the H2 molecule from a 
self- consistent solution of the OAO scheme explained in the 
text and compared to a hydrogenic (Is) orbital 4>h{z). Both 
curves are plotted along the electronic z axis. In atomic units. 



Fig. Elvisualizes the optimized atomic orbital as obtained 
from Eq. (|119|1 for the H2 molecular ion. Compared to 
a hydrogenic Is orbital, which is added to the plot in 
dashed linestyle, we clearly see the anticipated influence 
of the other nucleus: At a distance of i? = (i?) ~ 2.2 a.u., 
a second peak appears, leading to what we called a po- 
larized orbital. We may view this orbital as being opti- 
mized in the sense that it provides the best ground-state 
solution when used in the expression H124|) for S^". In- 
deed the results obtained for the H2 and H2 molecules 
molecule, which are again given in Tables HI and Hll consis- 
tently improve upon the Hartree data. This remains true 
for the harmonic constant w, where the deviations found 
in the Hartree scheme are somewhat reduced within the 
current approach. Correspondingly, the nuclear densities 
and potentials are slightly improved, as seen in Figs. 
-0] However, the disagreement with the exact curves is 
still quite large. 

To further investigate this point, we have, following 
Eq. (|125|l . decomposed the conditional potential into its 
different parts. The results obtained for the H2 molecule 
are shown in Fig.[7| where V^ond,w(^), Eq. H126|) . is plot- 
ted on the left-hand side and V^ond,T(^), Eq. (|127|l . is 
plotted on the right-hand side. In addition to the curves 
obtained from the OAO and the BO approach, which 
again serves as a reference, we added the results cal- 
culated from a simple LCAO ansatz using hydrogenic 
(Is) atomic orbitals in Ea. H114|l instead of the optimized 
atomic orbitals. The corresponding curve is denoted by 
H(ls)-LCAO in Fig.0 We first observe that the results 
from the simple LCAO and the optimized OAO scheme 
are very similar in shape. Yet, the ones obtained from the 
optimized OAO are very close to the exact (BO) numbers 
at the equilibrium internuclear distance {R), as clearly 
visible in Fig.[7| This enables the OAO approach to pre- 
dict ground-state properties nicely, whereas the simple 
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FIG. 7; Contributions to the nuclear conditional potential 
fl25}} as obtained from a self-consistent solution of the MCKS- 
OAO scheme for the molecule. They are compared to the 
corresponding BO curves and to results provided by a simple 
LCAO employing hydrogenic (Is) orbitals. In addition, the 
mean (equilibrium) internuclear distance {R) is marked. In 
atomic units. 



LCAO only leads to qualitatively correct results. How- 
ever, since the optimized OAO potentials are basically 
shifted (see FigEJ, the i? — > cx) asymptotics, which - by 
construction - is correctly described by the simple-OAO 
approach, is incorrect in the present approach. In the 
unified atom limit, on the other hand, both OAO schemes 
produce large errors. This is an inherent shortcoming of 
the OAO approach, which is not set up to satisfy the 
R^O limit. 

In conclusion, the results of the method presented in 
this section are clearly superior to the ones obtained from 
the Hartree approach. Since the atomic orbitals are opti- 
mized for the molecular ground state, the scheme works 
nicely for quantities depending mostly on the equilib- 
rium geometry but fails to substantially improve on the 
Hartree scheme for large internuclear distances. This de- 
ficiency will be dealt with in the next section. 



D. Scaling the Atomic Orbitals 

Owing to the successes of the OAO method to describe 
the bonding region, an ansatz similar to H114|) will lay the 
foundation also for the work presented in this section. 
However, in order to improve on the shortcomings of the 
OAO approach, we additionally concentrate on the task 
of setting up the scheme such that the separated as well 
as the unified atom limit are correctly reproduced. 

The following investigations are again based on the 
quantity which is considered to be the key quantity to 
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approximate, namely the electronic conditional density 
7(r|R). As above, we start out by decomposing the con- 
ditional density into orbital contributions as in Eg . 111 131) 
which are approximated by an LCAO-type ansatz 



7,(r|R) 



1 



and 



Ma /I 



(130) 



M„ 



Re, as before. The denominator 



Vj (R) arises from the normalization constraint (|115|l and 
is given by 

1 



dr 



(131) 



As a consequence, the first sumrule (|65() for the condi- 
tional density is automatically satisfied for any choice of 
the atomic orbitals. 

Up to this point, we just repeated the analysis of 
Sec. II V CI Now, we have to specify the atomic or- 
bitals {4>f^^}- In the approach presented in the last 
section, they were determined by optimizing the molecu- 
lar ground-state configuration. Since these orbitals were 
further used to describe the large-i? behavior of a di- 
atomic molecule, the deviations reported on above were 
found. 

In order to improve on that, we consider the separated 
atom [R oo) limit. There, the system consists of two 
atoms A and i3, which do not interact among each other. 
The atoms can thus be described by electronic orbitals, 
denoted by {4>J^j }, which yield the ground-state densi- 
ties Pa/b of the fragments. We assume that these orbitals 
are known, e.g., from an electronic DFT calculation for 
the single atoms. If the orbitals {4>^^} were inserted 
in Eq. (|130|l . we would obtain the simple LCAO scheme 
which was used for comparison in the last section. Obvi- 
ously, the bonding effects are not satisfactorily described 
within such a simple ansatz. In view of the mutual influ- 
ence of the atoms, we expect the orbitals to change when 
the atoms approach each other as illustrated by Fig.|B|in 
the last section. 

Here, we account for this change by introducing con- 
tracted orbitals 

0f ^(r) ^ 4'f{v) = Af ^^^/f (A,r). (132) 

The idea to model bonding effects by such a scaling pro- 
cedure can be explained in terms of the virial theorem: 
A decrease of the total energy due to chemical bonding 
leads to an increase of the kinetic energy and thus to 
a contraction of the orbitals js^]. Evidently, the size of 
the contraction depends on the molecular configuration 
and, in particular, on the internuclear distance. It should 
therefore be determined self-consistently from the densi- 
ties which characterize the system, i.e.. 



A, -A,[r,p](R) 



(133) 



At this point, we already see some benefits of this ap- 
proach. If, for the equilibrium geometry, a scaling param- 
eter A > 1 is used, the description of molecular bonding 



is improved upon the simple LCAO approach. On the 
other hand, employing A — > 1 for R — > oo, the large- 
R asymptotics of the conditional density is correctly re- 
produced. In the following, we will describe a way to 
calculate the scaling function Aj(R) within the MCKS 
scheme. We start by employing the second sumrule (|66ll . 
To construct our functional it is assumed that this equa- 
tion also holds true for its analogue formulated in terms 
of the orbital densities (as it does in the BO limit). 



dRr(R)7j(r|R), 



(134) 



where Pj{v) — |c/?j(r)p. If Eq. H134|) is satisfied for all j, 
the sumrule (|66|) is obeyed, too. Employing Eqs. H13U|) 
and H132f) . Eq. (|134|l is rewritten as 



p,(r)^/.Rr(R)|^ 

(/.^^,.(A,(R)r^)+<,.(A,(R)rB 



(135) 



Once the atomic orbitals {4''tl^} ^'"S given, the above in- 
tegral equation determines Aj(R) as an (implicit) func- 
tional of the nuclear density V and of the electronic 
MCKS orbital densities pj. However, a full solution of 
this integral equation is rather complicated and will not 
be attempted in the present approach. 

Instead, a simplified scheme appears highly desirable. 
To this end, we investigate the limits of the scaling func- 
tion Aj(R). As already noted above, we impose the con- 
dition that 



A,(R) 



(136) 



which guarantees the correct R ^ oo asymptotic behav- 
ior of the conditional density 7(r|R). Next, we consider 
the unified atom limit i? = 0. In that case, the orbital 
conditional density reads 



7,(r|R-0) 



Ag(0) 



0l^.(A,(O)r)+0l,(A,(O)r) 



(137) 

This quantity should be equal to the electronic orbital 
density of the unified atom, p^+Bj(r). We therefore 
choose Aj(0) such that 7(r|0) most closely resembles 
PA+s,i(i")- In other words, Aj(0) is obtained from the 
minimization problem 



mm 

A,(0) 



dr 



7i(r|0) -pA-^i3j(r) 



(138) 



where the unified-atom density is assumed to be known. 

We illustrate this prescription for H2 molecule: In the 
separated atom limit, the electron is represented by a 
hydrogenic (Is) orbital 



1 



: exp (- 



(139) 
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whereas the density of the unified - in this example: He~^ 
- atom reads 



PA+Air) = PHe+{r) 



exp{-2ZHe+r). (140) 



From Eq. H138|) . we immediately obtain 
A(0) = Z^,+ . 



(141) 



Therefore, A(0) is given by the sum of the charges of the 
two nuclei. In a more complicated system, we expect the 
bare nuclear charge Za+b to be replaced by an effective 
one. Employing Eq. (|141|l . it is easily seen that the condi- 
tional density (|130|l reproduces the correct unified atom 
limit. 

Furthermore, the small-i? behavior of the conditional 
density is analyzed. This can be done by expanding the 
electronic Hamiltonian He in powers of R, yielding 



He — Ha^ 



M1Z2-M2Z1 

Mn„r. 



^^i? + 0(i?2), (142) 



3 ■' 



where Ha+b denotes the Hamiltonian of the unified 
atom. For homonuclear systems, to which we restrict 
ourselves in all numerical calculations, the first-order cor- 
rection in Eq. H142(l vanishes. From the fact that the elec- 
tronic density corresponding to (|142ll basically coincides 
with 7(r|R), we obtain the small-i? {R 0) behavior: 



7(r|R)-p^+^(r)+0(i?2). 



(143) 



In view of the limits discussed above, we propose a 
simple parameterization of the scaling function: 



A,(R) = 1 



l + /3ji?^ 



(144) 



Using such a form, the R ^ 00 limit (|136|l is fulfilled. 
The constant aj follows from the unified atom limit, 
Eq. p38(l : aj — Xj{0) — 1. The exponent 7 is chosen 
such that the model conditional density behaves - for 
homonuclear molecules - as (|143|) for R ^ 0, leading to 
7 = 2. We are therefore left with one still unknown coef- 
ficient, namely Pj. To determine this constant, we resort 
to the integral equation (|134|l . Employing additionally 
quasi-classical nuclei, r(R) = (5(R — (R)), we obtain 



p,(r)=7,(r|(R)). 



(145) 



The coefficient Pj is then obtained self-consistently from 
fitting the model conditional density to the electronic 
MCKS orbitals densities: 



mm 

ft- 



dr 



p,(r)-7,(r|(R)) 



(146) 



Having calculated f3j, we put together all ingredients 
for the construction of the model conditional density. 



which is denoted by SAO (scaled atomic orbital) in the 
following, and finally arrive at 



y«^0[r,p](r|R) = 5:.^^'^'') 



2i.,(R) 



(147) 



0^^,.(A,(R)r^) +0^^^.(A,(R)rB 



with Aj (r) given by Eq. H144|l . Summarizing the above 
prescription, the parameters in Xj are obtained from (i) 
the atomic orbitals corresponding to the unified and the 
separate atom limit, which have to be provided as an 
input, and (ii) self-consistently via Eq. H14t)|) from the 
MCKS (orbital) densities. Thereby, the conditional den- 
sity 7^^'-'[r,p] is an (implicit) functional of the MCKS 
densities. Moreover, the SAO conditional density now 
satisfies, by construction, both normalization sumrules 
(the second sumrule at least in a good approximation), 
reproduces the correct asymptotic behavior and hence 
meets all the requirements set up in the beginning of the 
section. 

Having obtained an approximation for the conditional 
density, we again have to face the problem of the 
coupling-constant integration in Eq. Ht)2|) . One possibil- 
ity to overcome this problem was discussed in Sec. IIVCI 
where the He energy functional is expressed exclusively in 
terms of the conditional density at full coupling strength. 
Employing this expression, Eq. (|59|1 . we finally obtain 

i?ic°[r,p] = I dRr(R) (148) 

J dr WeniR,r) 7"^°(r|R) -f F^[7^^0](R) - 

The corresponding nuclear conditional potential is given 
by 

+ K'o^d°H(R) + K'on?xc(R). (149) 

with 

Kon°w(R) = y'driyen(R,r)7'^°(r|R) (150) 
t/cond,T(K) - 8 2^7 7fAO(r|R) 



- Ts.eip] (151) 



2 . _ 

- Elilp] (152) 
K'o^d°xc(R) = i?^c[7'^°](R)-i?xc[p]- (153) 

As above, Eq. (|149|l can be used in the MCKS scheme. 

As an interesting aside, we add an alternative ap- 
proach to calculate the nuclear conditional potential di- 
rectly from 7^=^(r|R). The idea rests on the observation 
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that T4ond(R-) is practically identical to the lowest energy 
BO-PES, if non-BO effects are negligible. Employing the 
Hellmann-Feynman (electrostatic) theorem we ob- 
tain 

^7(r|R) 



ayco„d(R) 



(9R 



= / dr 



(154) 



Evidently the slope of the nuclear conditional potential 
is solely determined by the conditional density at full 
coupling strength and could therefore be evaluated using 
Eq. H147|l . Moreover, compared to Eq. 1149(1 . the expres- 
sion l|154|) seems to be more efficient from a numerical 
point of view. However, the first approach leading to 
Eq. ((149|l proved to be more accurate and will therefore 
be used in the calculations presented below. 

To summarize, the nuclear conditional potential is cal- 
culated from Eq. (|149|l by using the model SAO con- 
ditional density (|147|l . As for the OAO approach, we 
furthermore neglect the nuclear response potential, ap- 
proximate the electronic potential by the Hartree ansatz, 
and solve the MCKS scheme self-consistently. 
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FIG. 9: Comparison of different contributions to the nu- 
clear potential, Eqs. HSUtl and H51\) . as obtained for the 
molecule from the BO, the SAO, and the OAO method. 
In atomic units. 



1 . Results 

First of all, we consider the scaling function \{R) ob- 
tained from a self-consistent MCKS calculation for the 
H2 molecule, which is plotted in Fig. |S1 By construc- 
tion, \{R) tends to one for large R such that the dis- 
sociation limit is correctly reproduced. For small R, on 
the other hand, we find that A(0) « 1.7. As expected 
above, this number is similar to the effective charge one 
obtains for the He atom within a simple Hartree-Fock 
treatment employing hydrogenic orbitals. At the equilib- 
rium distance (i?), the scaling function acquires a value 



1.6 



1.2 - 

X{<n>) 




1 <R> 2 3 4 5 6 7 

R 

FIG. 8: Scaling function \{R), Eq. \144\l , obtained from 
a self- consistent solution of the MCKS scheme for the H2 
molecule. In atomic units. 



of about 1.17, leading to a contraction of the orbital by 
this factor. As a consequence, the bonding energy is low- 
ered compared to a simple LCAO ansatz. Indeed, the 
ground-state energy of the H2 molecule obtained from 
the SAO approximation is close to the exact one, as is 
seen from the last column of Tab. The energy im- 
proved on the Hartree and on the LCAO data, reducing 
the error to about 0.5%. We also observe that the R- 
expectation values are slightly overestimated, which can 
be viewed as a left-over from the simple LCAO method, 
which generally tends to overestimate the bonding dis- 
tances. The remaining deviations can be attributed to 
changes in the orbital like, e.g., the appearance of a sec- 
ond peak as seen in Fig|^ which cannot be accounted for 
by the simple scaling procedure used in the SAO scheme 
presented here. This effect seems to be more pronounced 
for the iJj'^ molecule. From Tab.HTlwe find that especially 
the ground-state energy is somewhat worse than the re- 
sults obtained from the other approximations. However, 
we find a remarkable improvement in the harmonic con- 
stant. For both the H2 as well as the H2 molecule, the 
relative error in lo is lowered by more than an order of 
magnitude to about 3%. Correspondingly, the nuclear 
densities and potentials are expected to be closer to the 
exact results, too. From Figs. ^ and |3| we indeed find 
that the nuclear potentials obtained from the SAO ap- 
proach arc almost indistinguishable from the BO-PES 
in the asymptotic {R and R —^ 00) regime. Of 
course, this is hardly surprising, since the correct asymp- 
totic behavior was imposed in the construction of the 
conditional density. However, this consequently leads to 
a much improved shape of the nuclear potential, which 
is obvious from a comparison of the SAO curves to the 
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Hartree or LCAO data and is also reflected in the better 
harmonic constant reported above. Moreover, consider- 
ing the H2 molecule, the SAO nuclear potential coin- 
cides with the BO-PES not only asymptotically, but in 
the whole R range. This then leads to a nuclear density, 
plotted in Fig. 01 which nicely agrees with the BO one. 
For the H2 molecule, the agreement at intermediate in- 
ternuclear distances is not as good. From Fig.|21 we find 
that the nuclear density is slightly shifted to larger R. 
Additionally, the minimum of the SAO nuclear potential 
is too high, leading to the deviations in the energy al- 
ready mentioned in the discussion of Tab. ^ Still, the 
SAO approach provides the best overall description also 
for the H2 molecule. In particular, the nuclear density 
obtained from the SAO approach is closest to the exact 
one. Furthermore, as seen from Fig |^ where the dif- 
ferent contributions to the nuclear conditional potential, 
Eqs. (|15()|l and H151|l . are shown, the SAO curves are in 
satisfactory agreement with the BO results in the whole 
range of internuclear distances and thus improve on the 
OAO method, which only reproduces the correct values 
around the equilibrium distance. We therefore observe 
the effect of incorporating the unified as well as sepa- 
rated atom limit into the construction of the He energy 
functional. At this point, we also emphasize the impor- 
tance to account for the additional contributions which 
arise from the coupling-constant integration. As is seen 
from the right-hand side of Fig the i?-dependence of 
these terms is significant, and it would not be a good 
approximation to replace by Wen- In conclusion we 
find that the SAO approximation to the electron-nuclear 
correlation functional gives a very good overall descrip- 
tion of the BO potential. 

V. CONCLUSIONS 

For a unified quantum mechanical treatment of nuclear 
and electronic degrees of freedom we extended the tradi- 
tional density functional method to multicomponent sys- 
tems. We first discussed the choice of appropriate densi- 
ties serving as fundamental variables of the theory. It was 



shown that the usual definition of single-particle densities 
in terms of inertial coordinates is not well suited for the 
purpose of this work because such densities, as a conse- 
quence of Galilean invariance, are constant for all isolated 
systems and therefore not characteristic for their inter- 
nal properties. A suitable set of densities was obtained 
by defining the electronic density with respect to a coor- 
dinate system attached to the nuclear framework whereas 
the nuclear degrees of freedom were described by the di- 
agonal of the nuclear A^„-body density matrix. For these 
fundamental variables the Hohenberg-Kohn theorem and 
the Kohn-Sham equations were derived and the corre- 
sponding density functionals were analyzed in detail. The 
main new ingredient of the multicomponent theory is 
the electron- nuclear correlation functional. For this func- 
tional several approximations were derived and tested on 
the H2 and H2 molecules. It was found that the simplest 
Hartree approximation fails to give a good description of 
the bonding curve of these molecules. Considerable im- 
provement was obtained using an approximation based 
on optimized atomic orbitals. This method still had some 
deficiencies in the large and small bond distance limits. 
These deficiencies were finally removed using an approx- 
imation based on scaled atomic orbitals. Based on this 
first experience with MCDFT we can say that it presents 
a promising new approach to study electron-nuclear cor- 
relation phenomena beyond the BO approximation. A 
promising new field of applications seems the first prin- 
ciple treatment of electron-phonon interactions within a 
linear response language using a time-dependent exten- 
sion of the present theory. Another field of future ap- 
plications will be the study of combined ionization and 
dissociation dynamics of molecules in strong laser fields. 
For this case some approximations of similar spirit as 
discussed in this paper have already been applied succes- 
fully [HEI. 

VI. ACKNOWLEDGEMENTS 

RvL acknowledges support of 'Stichting voor Funda- 
menteel Onderzoek der Materie (FOM)'. 



[1] P.Hohenberg and W.Kohn Phys.Rev.136, B864, (1964) 
[2] W.Kohn and L.J. Sham Phys.Rev.140, A1133, (1965) 
[3] T. Kreibich and E.K.U. Gross. Phys. Rev. Lett., 86, 
2984, (2001). 

[4] R.van Leeuwen Phys.Rev. B69, 115110 (2004); 
199901(E) (2004). 

[5] M. Liiders, M.A.L. Marques, N.N. Lathiotakis, A. Floris, 
G. Profeta, L. Fast, A. Continenza, S. Massidda, E.K.U. 
Gross, Phys. Rev. B 72, 024545 (2005). 

[6] M.A.L. Marques, M. Liiders, N.N. Lathiotakis, G. Pro- 
feta, A. Floris, L. Fast, A. Continenza, E.K.U. Gross, S. 
Massidda, Phys. Rev. B 72, 024546 (2005). 

[7] A. Floris, G. Profeta, N.N. Lathiotakis, M. Liiders, 



M.A.L. Marques, C. Franchini, E.K.U. Gross, A. Con- 
tinenza, and S. Massidda, Phys. Rev. Lett. 94, 037004 

(2005) . 

[8] G. Profeta, C. Franchini, N.N. Lathiotakis, A. Floris, A. 
Sanna, M.A.L. Marques, M. Liiders, S. Massidda, E.K.U. 
Gross, A. Continenza, Phys. Rev. Lett. 96, 047003 

(2006) . 

[9] G.D.Mahan Many-Particle Physics Plenum Press (NY), 
(1990). 

[10] K.Hannewald and P.A.Bobbert Phys.Rev. B69, 075212 
(2004). 

[11] Z.Bacic and J.C.Light Ann.Rev.Phys.Chem. 40, 469 
(1989). 



23 



[12] D.Dulic, S.J.van der Molen, T.Kudernac, H.T.Jonkman, 
J.J.D.de Jong, T.N.Bowden, J.van Esch, B.L.Feringa and 
B.J.van Wees Phys.Rev.Lett. 91, 207402 (2003). 

[13] T.Kreibich Multicomponent Density-Functional Theory 
for Molecules tn Strong Laser Fields, Ph.D. Thesis, Uni- 
versity of Wurzburg (2000) 

[14] T.Kreibich, N.I.Gidopoulos, R.van Leeuwen and 
E.K.U. Gross in Progress in Theoretical Chemistry 
and Physics, VoL14, "The Fundamentals of Electron 
Density, Density Matrbc and Density Functional The- 
ory in Atoms, Molecules and the Solid State" Eds. 
N.I.Gidopoulos and S.Wilson (Kluwer, Dordrecht) 2003. 

[15] T.Kreibich, R.van Leeuwen and E.K.U. Gross, 
Chem.Phys. 304, 183 (2004). 

[16] J.F. Capitani, R.F. Nalewajski, and R.G. Parr. 
J. Chem. Phys., 76 568, (1982). 

[17] H. Goldstein. Classical Mechanics. Addison- Wesley, 
Reading, 1980. 

[18] F.Villars Nucl.Phys. 3, 240 (1957) 

[19] F.M.H.Villars and G.Cooper Ann.Phys. 56, 224 (1970) 
[20] B.Buck, L.C.Biedenharn, and R.Y.Cusson Nucl.Phys. A 

317, 205 (1979) 
[21] C.Eckart Phys.Rev. 47, 552 (1935) 

[22] E.B.Wilson, J.C.Decius and P.C.Cross, Molecular Vibra- 
tions (Dover, New York, 1955) 

[23] J.D.Louck and H.W.Galbraith Rev.Mod.Phys.48, 69 
(1976) 

[24] J.D.Louck J.Mol.Spectrosc. 61, 107 (1976) 

[25] P.R. Bunker. Molecular Symmetry and Spectroscopy. 

Academic Press, New York, 1979. 
[26] B.Sutcliffe J.Chem.Soc, Faraday Trans. 89, 2321 (1993) 
[27] B.Sutcliffe Adv. Chem.Phys. 114, 1 (2000) 
[28] R.G.Littlejohn and M.Reinsch Rev.Mod.Phys. 69, 213 

(1997) 

[29] R. Schinke. Photodissociation Dynamics. Cambridge 

University Press, Cambridge, 1993. 
[30] L.M. Sander, H.B. Shore, and L.J. Sham. 

Phys. Rev. Lett. 31 533, (1973). 
[31] R.K. Kalia and P. Vashishta. Phys. Rev. B17, 2655, 

(1978). 

[32] N. Gidopoulos. Phys. Rev. B57, 2146, (1998). 
[33] X.Lopez, J.M.Ugalde and E.V.Ludeiia Eur.Phys.J. D37, 
351 (2006) 

[34] R.M. Dreizler and E.K.U. Gross. Density Functional 

Theory. Springer, Berlin, 1990. 
[35] M. Levy Proc. Natl. Acad. Set. USA 76, 6062, (1979). 
[36] G.Chaban, J.O.Jung and R.B.Gerber J. Chem.Phys. Ill, 



1823 (1999) 

[37] J. Harris and R.O. Jones. J. Phys. F4, 1170, (1974). 
[38] D.C. Langreth and J. P. Perdew. Solid State Commun. 

17 1425, (1975). 
[39] O. Gunnarsson and B.I. Lundqvist. Phys. Rev. B13, 

4274, (1976). 

[40] G. Hunter Int. J. Quant. Chem. Symp. 8, 413 (1974) 

[41] G. Hunter Int. J. Quant. Chem. 9, 237 (1975) 

[42] G. Hunter Int.J.Quant.Chem.l7, 133 (1980) 

[43] G. Hunter Int.J.Quant.Chem.l9, 755 (1981) 

[44] J.Czub and L.Wolniewicz Mol.Phys. 36, 1301 (1978) 

[45] P.Cassam-Chenai Chem.Phys.Lett. 420, 354 (2006_2_ 

[46] N.I.Gidopoulos and E.K.U. Gross cond-mat/0 502433| 
(2005) 

[47] M. Berry. Proc. Phys. Soc. (London), A392, 45 (1984). 
[48] C. Lanczos. J. Res. Nat. Bur. Stand. 45, 255, (1950). 
[49] I.L. Thomas. Phys. Rev. 185, 185, (1969). 
[50] I.L. Thomas and H.W. Joy Phys. Rev. A2, 1200, (1970). 
[51] R.T. Sharp and G.K. Horton. Phys. Rev. 90 317, (1953). 
[52] J.D. Talman and W.F. Shadwick. Phys. Rev. A14, 36, 
(1976). 

[53] T. Grabo, T. Kreibich, S. Kurth, and E.K.U. Gross. 
Strong Coulomb Correlations in Electronic Structure: 
Beyond the Local Density Approximation, pages 203-311. 
Gordon&Breach, Tokyo, 1999. 

[54] R. van Leeuwen. PhD thesis, Vrije Universiteit, Amster- 
dam, 1994. 

[55] R.P. Feynman. Phys. Rev. 56, 340, (1939). 

[56] E.H. Lieb. In A. Shimony and H. Feshbach, editors, 
Phystcs as Natural Philosophy, page 111. MIT Press, 
Cambridge, 1982. 

[57] A. Shapere and F. Wilczek. Geometric Phases in Physics. 
World Scientific, Singapore, 1989. 

[58] R.van Leeuwen and E.K.U. Gross in "Time-Dependent 
Density Functional Theory" Eds. M. A. L. Marques, 
C.A.Ullrich, F.Nogueira, A.Rubio, K.Burke and 
E.K.U. Gross, Lectures Notes in Physics Vol.706, 
Springer (Berlin) (2006) 

[59] In standard electronic DFT, one can prove that the min- 
imum of F exists IH^ . 

[60] In the formalism presented here, the Berry-phase effects 
are assumed to be representable by a scalar nuclear po- 
tential. In view of the connection of the Berry phase to 
a vector potential in the nuclear Schrodinger equation 
fH^ . a multicomponent current- density description ap- 
pears more appropriate to treat these effects. 



